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OJQi ABSTRACT 

Context. Theoretical spectra of X-ray bursting neutron star (NS) model atmospheres are widely used to determine the basic NS pa- 
rameters such as their masses and radii. Compton scattering, which plays an important role in spectra formation at high luminosities, 
is often accounted for using the differential Kompaneets operator, while in other models a more general, integral operator for the 
Compton scattering kernel is used. 

Aims. We construct accurate NS atmosphere models using for the first time an exact treatment of Compton scattering via the integral 
relativistic kinetic equation. We also test various approximations to the Compton scattering redistribution function and compare the 
results with the previous calculations based on the Kompaneets operator. 

Methods. We solve the radiation transfer equation together with the hydrostatic equilibrium equation accounting exactly for the ra- 
i-C 1 diation pressure by electron scattering. We use the exact relativistic angle-dependent redistribution function as well as its simple 

approximate representations. 

Results. We thus construct a new set of plane-parallel atmosphere models in local thermodynamic equilibrium (LTE) for hot NSs. 
The models were computed for six chemical compositions (pure H, pure He, solar H/He mix with various heavy elements abundances 
1 Z = 1, 0.3, 0.1, and 0.01Z Q , and three surface gravities logg = 14.0, 14.3, and 14.6. For each chemical composition and surface 

gravity, we compute more than 26 model atmospheres with various luminosities relative to the Eddington luminosity L Edd computed 
for the Thomson cross-section. The maximum relative luminosities L/Leu reach values of up to 1.1 for high gravity models. The 
emergent spectra of all models are redshifted and fitted by diluted blackbody spectra in the 3-20 keV energy range appropriate for the 
RXTE/PCA. We also compute the color correction factors f c . 
Conclusions. The radiative acceleration g rad in our luminous, hot-atmosphere models is significantly smaller than in corresponding 
models based on the Kompaneets operator, because of the Klein-Nishina reduction of the electron scattering cross-section, and there- 
\Q . fore formally "super-Eddington" model atmospheres do exist. The differences between the new and old model atmospheres are small 

for L/Leu < 0.8. For the same g m d/g, the new f c are slightly larger (by approximately 1%) than the old values. We also find that the 
model atmospheres, the emergent spectra, and the color correction factor computed using angle-averaged and approximate Compton 
scattering kernels differ from the exact solutions by less than 2%. 

Key words, radiative transfer - scattering - methods: numerical - stars: atmospheres - stars: neutron - X-rays: stars 



1. Introduction body dGallowav et al.l I2008I). Theoretica l models of hot NS 



X 



atmospheres show (London et al.l Il986t, lLapidus et al.1 Il986: 



. X-ray bursting neutron stars (NSs) are members of low- Pavl ov et all 1 19911) that the emergent spectra are close to di- 
. mass X-ray binaries with quasi-periodi cal thermonuclear mte( i blackbody spectra F E * f<: A B E (f c T eSS ) owing to strong 
flashes on their surfaces (s ee reviews by iLewin et all 119931: energy exchange caused by Compton scattering. The color cor- 
IStrohmaver & Bildstedl2006l). Thermonuclear burning occurs at rection factor f c = T Q /T eS , which relates the color tempera- 
the bottom of the freshly accreted matter and can be so pow- ture of the spectrum T c to the effective temperature of the at- 
erful that the luminosity reaches the Eddington limit. These mosphere r eff , takes values of about 1.3-1.9. Theoretical de- 
bursts lead to photospheric radius expansion (PRE) and are a pendences of the color correction factor on luminosity for var- 
potent ially powerful tool for determining t h e NS masses and ious chemical compositio ns and gravities were computed by 
radii dEbisuzakil 119871: iDamen et all 119901: Ivan Paradiis et all ISuleimanov et al.l (l2011bl hereafter SPW11). The atmosphere 
119901). The knowledge of NS basic parameters is extremely mo dels used in their wo rk were based on the Kompaneets opera- 
important for establishing the physical properties (equation of tor dKompaneetsll 19571) description of Compton s cattering. This 



state) of supra-nuclear de nse matter in the NS inner cores (see approach was also used in most previous s t udies dLondon et al.l 

lLattimer & Prakash||2007|, for a review). 19861 lLapidus et al.lll986t lEbisuzakilfl987t iPavlov et al.lll991l) . 

The precise radius determination of the NSs from X-ray The Kompaneets operator describes Compton scattering in 

bursts is impossible without accurate spectral models. The ob- the non-relativistic, isotropic, diffusion approximation, which 

served spectra of X-ray bursts are often well-fitted by a black- should still be rather adequate for hot NS model atmospheres 
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with typical effective temperat ures below ~2 ke V. On the other 
hand, Madej and coll aborators (lMadeilll99ll iMadei et al1l2004t 
Maic zvna et al.ll2005l) used the integral description of Compton 
scattering empl oying the angle- averaged redistribution function 
(RF) derived bv lGuilbertl Jl98lb . A fully relativistic treatment of 
Compton scattering becomes important at luminosities close to 
the Eddington values, because of the reduction of the effective 
cross-section and the corresponding drop in the radiation pres- 
sure force. This obviously changes the maximum value of the lu- 
minosity where hydrostatic equilibrium can still be achieved, i.e. 
the actual value of the Eddington luminosity. This has important 
implications for the method of determining the NS masses and 
radii from PRE bursts, which is bas ed on the Eddington lumi- 
nosity (see Suleim anov et al.|[201 lal hereafter Sll). It is there- 
fore, necessary to check the validity of model atmospheres and 
spectral properties computed with the Kompaneets operator by 
comparing the results with more accurate models based on an 
exact treatment of Compton scattering using relativistic kinetic 
equations. 

Here we construct new models of hot NS atmospheres based 
on an exact treatment of Compton scatte ring that employs a 
fully re lativistic, an gl e-depe nd ent RF jA haronian & Atovan 
198 U iPrasad et alj 1 19861: iNagirner & Poutanenl [1994; 
Poutanen & Svensson 1996). In Sect. [2] we present the 



main equations describing the atmosphere models as well as 
the methods for their solution. We also compare atmosphere 
models based on various RFs for Compton scattering and using 
the Kompaneets operator. In Sect. [3] we present a new set of 
atmosphere models, emergent spectra, and the color correction 
factors. We compare them to previous results obtained with 
the Kompaneets operator. In Sect. |4] we discuss applications 
of the new models to the observations of X-ray bursts and 
determination of the NS masses and radii. We summarize in 
Sect. [5] In Appendix [A] we derive the exact RF for Compton 
scattering. The details of the method of solution of the radiative 
transfer equation are presented in Appendix [B] Comparison 
with the previous attempts to compute atmosphere models 
using integral approach to Compton scattering are given in 
Appendix [C] And, finally, Appendix [D] presents the spectral 
characteristics of the models. 

2. Method of atmosphere modeling 

2.1. Main equations 

A method for modeling hot X-ray bursting NS atmospheres that 
employs the Kompaneets operator for Compton scattering was 
described in detail in SPW1 1. Here we repeat the basic assump- 
tions and emphasize the differences arising because of the new 
treatment of Compton scattering as well as radiation transfer. 

We compute models in hydrostatic and radiative equilibrium 
in a plane-parallel approximation. The main input parameters 
are the chemical composition (particularly the hydrogen mass 
fraction X), the surface gravity 

GM 

and the relative NS luminosity / = L/Lndd, where LEdd is the 
Eddington luminosity measured at the NS surface 



\tiGMc 

^Edd = (1 +Z), 



for the Thomson scattering opacity 



2 „-l 



* e = cr T — « 0.2(1 + X)cm / g 
P 



(2) 



(3) 



Here <x T = 6.65 x 10~ 25 cm 2 is the Thomson cross-section, p 
is the gas density, and N e is the electron number density. The 
gravitational redshift is related to the NS parameters as 



1 +z = (1 -2GM/c 2 R) 



-1/2 



(4) 



We assume that the flux is constant throughout the NS surface. 
Our calculations are valid for a patch on the NS surface if in- 
stead of the relative luminosity we consider the relative flux, 
I = F/F^ dd , as a parameter, where 



L Edd GMc 



Edd 



(1+z). 



AnR 2 R 2 K e 

The effective temperature T e ff can be expressed via / as 



(5) 



(6) 



where the Eddington temperature Tndd is the maximum possible 
effective temperature on the NS surface, which is evaluated using 
the Thomson scattering opacity 



o~sbT, 



Edd 



•^Edd 



Ke ' 



(7) 



The structure of the atmosphere for an X-ray bursting NS is 
described by a set of differential equations. The first one is the 
hydrostatic equilibrium equation 



dm 



(8) 



where g ra d is the radiative acceleration and P g is the gas pressure 
and the column density m is defined as 



dm — —pds , 



(9) 



where s is the vertical distance. 

The second equation is the radiation transfer equation for the 
specific intensity I(x, fi) accounting for Compton scattering (see 
Appendix|A]for derivation). In the plane-parallel approximation, 
it has the form 



H-r- r = I{x,fi) - S (x,/j.), 

dT(X,/J.) 

where 

dr{x,jx) — [cr(x,/j) + k(x)] dm, 



(10) 



(11) 



p = cos 6 is the cosine of the angle between the surface nor- 
mal and the direction of radiation propagation, x = hv/m e c 2 
is the photon energy in units of electron rest mass, and k(x) is 
the "true" absorption opacity. The electron scattering opacity ac- 
counting for the induced scattering is 



cr(x,p)-K e — ^ x\dx\ j"dfi\R(xi,fi\;x, 



AO 



1 + 



ci(xi,m) 



and 



C 



1 

2m„ 



(12) 



(13) 



The source function is a sum of the thermal part and the scatter 
ing part 

k(x) 



S(x,fi) 



B, 



<r(x, jX) + k(x) * <t(x, jS) + k(x) 

dfj.iR(x,fi;xi,fii)I(xi,fJ.i) 



(14) 



Jo xi J-i 
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where 

dv 

B x = By- 
ax 



1 



C exp(;t;/0) - 1 
with By being the Planck function and 







kT 



(15) 



(16) 



the electron temperature. 

The RF R(x,/j,; xi,fi\) describes the probability that a photon 
with the dimensionless energy x\ propagating in the direction 
corresponding to [i\ is scattered to energy x and in a direction 
corresponding to yu. This function is found by integrating over 
the azimuthal angle <p of the RF R(x, xi,rf), which depends on the 
cosine of the angle between the directions of the photon propa- 
gation before and after scattering r\ 



R{x,n\x\,\x{) 



r*27T 



R(x, x\,T]) dip, 



(17) 



The RF depends on the depth s through the electron temperature 
and satisfies the relation (see Eq. IA.7t 



R{x\,ji\; x,ji) — R(x,n; xi,yui) exp 



(18) 



which is the consequence of the d etailed balance relation 
dPomraningl 1 9731: iNagirner & Poutanenl 19941 see AppendixlAt. 
This implies that the source function given by Eq. (TBI equals 
the Planck function for the photon field described by the Bose- 
Einstein distribution of any chemical potential. 

In this paper, we use two RFs for Compton scattering. The 
first one is the fully relativistic ex act RF valid for any photon 



1981 



y_p 

energy and electron temperature (Aharonian & Atovan 
Nagirner & Pou tane n 199 3, 1994; iPoutanen & Svenssonl I 
see Eqs. IA. 181 and IA.22I in Appendix IA.21 . The second one 



1996 



is an approximate RF corresponding to the isotropic scatter- 
ing in the electron rest frame (Eq. IA.26I in Appendix IA.21 . 
which is accurate at temperat ures below about 100 keV and non- 
relativistic pho t on energies (Arutyunyan & Nik ogosvan] Il980t 
lPoutanenlll994t|Poutanen & Svenssonlll996l) . In both cases, we 
also consider the angle-averaged RFs 



R{x,x\) 



I f +1 
2j-i 



drj R(x, x\,rj), 



substituting it instead of the angle-dependent RF into Eq. (fTTb . 

The formal solution of the radiation transfe r equati on (TTOb is 
obtain ed using the short-characteristic method (lOlson & Ku nasz 
fl987h in three angles in each hemisphere, and the full solution 
is found with an accelerated A-iteration method (see details in 
AppendixlBl. 

The radiation pressure acceleration g m & is computed using 
the RF as 



?rad = 



dPiad 27T d 

dm 



2n d r° r 1 . 

= r~ I ax i fi I(x, 

c dm Jo J_i 



fi)dfi 



(20) 



2tt C°° 

= — I dx I [<t(x, yu) + k(x)] [I(x, yu) - S (x, yu)] // d/j, 
c Jo J-i 

where the derivative with respect to m is replaced by the first 
moment of the radiation transfer equation ( TTOt . When the source 



functions and the opacities are isotropic, this expression is re- 
duced to the standard definition 



grad 

where 



_ 4n r 

c Jo 



[o~(x) + k(x)] H x (m) dx, 



yU I(X, yli) d\l 



(21) 



(22) 



is the first moment of specific intensity. 

These equations are completed by the energy balance equa- 
tion 



dx J [cr(x, //) + k(x)] [I(x, fx)-S (x, yu)] dfx = 0, (23) 
the ideal gas law 



P a = Mot kT, 



(24) 



where iV to t is the number density of all particles, and the parti- 
cle and charge conservation equations. In our calculations, we 
assumed local thermodynamic equilibrium (LTE), therefore the 
number densities of all ionization and excitation states of all ele- 
ments were calculated using the Boltzmann and Saha equations. 
We accounted for the pressure ionization effects on hydrogen 
and helium populations using t he occupation probabili ty formal- 
ism (|Hummer & Mihalasl[l988l) as described by Huben v"et al.l 
(11994 ). In addition to electron scattering, we took into account 
the free-free opacity as well as the bound-free transitions for all 
ions of the 15 most abundant chemical element s (H, He, C, N, 
O, Ne , Na, Mg, Al, Si, S, Ar, Ca, Fe, Ni) (see libra gimov et all 
120031) using opacities from lVerner & Yakovlevl ( 1995). 



2.2. Method of solution 

To solve the above equati ons, we used our version of the 
computer code ATLAS dKuruczl 1 19701 1 1993b. modified to 
deal with high temperature s dSuleimanov & Poutanenl [2006; 
ISuleim anov & Wernerll2007l) . The code was further developed 
to account for Compton scattering using the RF approach. 

In our computations, we used 300-360 logarithmically 
equidistant frequency points in the range 10 14 -10 20 Hz 4x 
10" 4 -400 keV) for the luminous model atmospheres (I > 0.1), 
and 10 14 — 10 19 Hz for I < 0.1. The calculations were performed 
at a set of 98 depth points m\ distributed equidistantly on the 
logarithmic scale from 10 -6 to m m!lx = 10 6 g cm~ 2 . The ap- 
nm propriate value of m max was chosen to satisfy the condition 
V T v,b-f,f-f('«max)Tv('«max) > 1 at all frequencies, where r Vjb -f,f-f 
is the optical depth computed with the true opacity only (bound- 
free and free-free transitions, without scattering). This require- 
ment was necessary for satisfying the inner boundary condition 
of the radiation transfer problem. 

The course of the calculations was the same as for the 
method that adopts the Kompaneets operator (SPW11). First, a 
starting grey atmosphere model was calculated and opacities at 
all depth points and all frequencies were obtained. The solution 
of the radiative transfer equation ( TTOb was checked for the energy 
balance equation ( f23l . together with the surface flux condition 



4n H x (m = Q)dx = AnH 
Jo 

The relative flux error 
H 



eft'- 



s H (m) = 1 



HJm)dx 



(25) 



(26) 
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Fig. 1. Emergent spectrum (top panel) and temperature struc- 
ture (bottom panel) of the fiducial model (pure hydrogen, T e g 
1 .8 X 10 7 K, logg = 14.0) computed using three different meth- 
ods. In accurate method 1 every A-iteration starts from the ther- 
mal part of the source function (solid curves give the results 
for the relative accuracy of 10" 4 ). In accelerated method 2, the 
A-iterations start from the source function taken from previous 
temperature iteration, but at every fifth temperature correction 
they start from the thermal part of the source function (dashed 
curves show results for the relative accuracy of 1CT 4 ). Method 
3 is the same as method 2, but for the relative accuracy of 10~ 5 
(dot-dashed curve). In the top lower panel, the ratios of the spec- 
tra for methods 2 and 3 to the spectrum computed with method 
1 are shown. The ratio of the temperature structures computed 
using methods 2 and 1 is shown in the bottom lower panel. 



and the energy balance error 



s A (m) 



= - I dx I [o-(x, p) + k(x)] [I(x, p)-S (x, p.)] dp (27) 
2 Jo J -i 



were calculated as functions of depth. Temperature corrections 
were then evaluated using three different procedures. In the up- 
per atmospheric layers, we used the integral A-iteration method, 
modified for Compton scattering, based on the energy balance 
equation ( f23l . The temperature correction for a particular depth 
was found as 



A7\ 



-e A (m) 



A d (x) - 1 



1 - a(x)Ad(x) 



k(x) -j^- dx 



(28) 



where a(x) = crcs(x)/(k(x) + crcs(x)), and Ad(jc) is the diagonal 
matrix element of the A-operator. Here cr cs (jc) is the Compton 
scattering opacity averaged over the r elativistic Maxwellian 
electron distribution (see Eq. (A16) in iPoutanen & Svenssonl 



119961 which is equivalent to Eq. (fT2l if one ignores the induced 
scattering). In the deep layers, we used the Avrett-Rrook flux 
correction based on the relative flux error sn(m). Finally, the 
third procedure w as the surface correction based on the emer- 
gent flux error (see iKuruczl 1 9701 for a detailed description of the 
methods). 

The iteration procedure is repeated until the relative flux er- 
ror is smaller than 0.1%, and the relative flux derivative error is 
smaller than 0.01%. As a result, we obtain a self-consistent NS 
model atmosphere, together with the emergent spectrum of ra- 
diation. We note that this accuracy is unachievable for luminous 
models with g ra d » g, and that these models can have larger rel- 
ative flux errors, up to 2-3%. 

2.3. Accuracy of computation 

To compute a new extended set of hot NS model atmospheres, 
we accelerated the convergence of the iterations of the radia- 
tion transfer equation by using the source function from the 
previous temperature iteration as the first approximation (see 
Appendix iBli. However, in every fifth temperature iteration the 
radiation transfer equation was solved using the pure thermal 
source function as a first approximation. We compared a pure 
hydrogen model atmosphere computed for r e ff = 1.8xl0 7 Kand 
logg = 14.0 (the fiducial model) using this accelerated approach 
with the model computed without acceleration. The temperature 
structures differ by less than 0.3 %, and the differences between 
the emergent spectra are about 1 % in the 3-20 ke V energy range 
(typical of RXTE/PCA) and larger in the Wien tail (Fig. 0Q>. 

As a convergence criterion for the solution of the radiation 
transfer equation, we chose the maximum relative error of 10~ 4 
in the mean intensity at all depths and energies. To determine 
the uncertainty in the final spectrum caused by this criterion, we 
compared the emergent spectra computed for the same model at- 
mosphere with the accuracies of 10~ 4 and 10~ 5 (see top panel of 
Fig.[TJ. We see that the relative error is smaller than 1 % at all en- 
ergies, which is then the intrinsic accuracy of our model spectra. 
We note that a similar error is introduced into the angular de- 
pen dence of the specific i ntensities by ignoring polarization (see 
e. p. IChandrase khar 1960, compare his Tables XV and XXIV). 

2.4. Various RFs and the Kompaneets operator 

A comparison of model atmospheres, which were computed us- 
ing our new code with four different RFs (exact angle-dependent, 
exact angle-averaged, approximate angle-dependent, and ap- 
proximate angle-averaged) is shown in the left panels of Fig. [2] 
The model computed using the approximate angle-dependent 
RF is almost indistinguishable from the reference model cal- 
culated using the exact angle-dependent RF. This is unsurpris- 
ing, because the approximate function matches the accurate one 
very well up to the tem peratures of about 10 keV Similar 
results were obtained by IPoutanen & Svenssonl d 19961) for the 
Comptonization spectra in the optically thin slabs. The differ- 
ences are smaller than 1 % for the emergent spectra and below 
0.3% for the temperature structure. Deviations of model atmo- 
spheres computed using the angle-averaged RF of the reference 
model are more significant, at about 2% for both the temperature 
structure and spectra. 

The model atmosphere calculated employing the 
Kompaneets operator has more significant differences from the 
reference model (see right panels in Fig. of up to 10% in 
the temperature structure and about 3^4% in the emergent flux 
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Column density (g cm" ) Column density (g cm" ) 



Fig. 2. Left panels present the emergent spectrum (top panel) and the temperature structure {bottom panel) of the fiducial model 
computed using four different RFs: (a) exact angle-dependent (solid curves), (b) exact angle-averaged (dashed curves), (c) approx- 
imate angle-dependent (dot-dashed curves), and (d) approximate angle-averaged (dotted curves). In the top lower sub-panel, the 
ratios of the spectra (b), (c), and (d) to the spectrum for case (a) are shown. In the bottom lower sub-panel, the ratio of the tem- 
perature structures (b), (c), and (d) to that of case (a) are shown. Right panels present the emergent spectrum (top panel) and the 
temperature structure (bottom panel) of the fiducial model computed using an exact angle-dependent RF (solid curves) and the 
Kompaneets operator (dashed curves). 



at 1 keV. The deviations are much smaller (< 1%) near the 
spectral maximum and larger in the Wien tail. A rather good 
agreement between model atmospheres computed with the 
relativistic exact angle-dependent RF and the non-relativistic 
angle-independent Kompaneets operator again is unsurprising, 
because temperatures of the upper atmosphere layers where 
the emergent spectra form are sufficiently low (~2-4 keV) and 
relativistic corrections are small. 



3. New grid of models 

3. 1 . General properties 

We computed a new set of hot NS model atmospheres using the 
exact relativistic angle-dependent RF. The models were calcu- 
lated for six chemical compositions (pure hydrogen, pure he- 
lium, and solar hydrogen/helium mix with various heavy ele- 
ment abundances: of solar and 0.3, 0.1 and 0.01 of solar). For 
every chemical composition, 26-28 models with relative lumi- 
nosities spanning the interval from / = 0.001 to 1.06-1.10 for 
three values of the surface gravity (\ogg =14.0, 14.3 and 14.6) 
were calculated. Because the radiative acceleration in our mod- 
els is smaller than in the models based on the Thomson opacity 
owing to the Klein-Nishina reduction in the cross-section (see 
below), there exist formally "super-Eddington" (relative to Ledd) 
models. The Klein-Nishina reduction depends on the electron 



temperature, which is higher for larger surface gravities, there- 
fore the limiting luminosity is higher for larger log g. 

Examples of emergent spectra and temperature structures for 
the models with logg = 14.0 and various chemical compositions 
(pure hydrogen and solar mix) are shown in Fig. [30 The cor- 
responding emergent spectra and temperature structures com- 
puted with the old code employing the Kompaneets operator are 
also shown. At low luminosities, the models are very close to 
each other, which is expected as at low temperatures the dif- 
fusion (Kompaneets) approximation is an accurate representa- 
tion of Compton scattering. For models close to the Eddington 
limit, the treatment of the radiative acceleration becomes im- 
portant. Contribution of the radiation force to the hydrostatic 
equilibrium g ra d/g is smaller in the new models despite I being 
the same for both model sets. It is well-known that the model 
spectra with large contributions of radiative ac celeration are 
harder, and their surface temperatures are higher ( London et alj 
119861: iLanidus et alJll986t lEbisuzakil[l987t iPavlov et alJl!991h . 
In complete agreement with this, the old models with large I are 
harder and hotter than new models. However, for the same g m d/g 
the new model spectra are hotter. Normalizing the spectra to the 
maximum flux and plotting them against the scaled photon en- 



1 The spectral energy distributions of fluxes and specific intensities 
for all models are described in Appendix iDl 
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Fig. 3. Comparison of the emergent spectra (top panel) and the temperature structures (bottom panel) computed using the present 
code with the exact RF (solid curves) and the previous code employing the Kompaneets operator (dashed curves). The models for 
various relative luminosities are marked by corresponding numbers next to the curves. Left panels are for pure H and right panels 
are for the solar abundances. 
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Fig. 4. Comparison of the emergent spectra (top panel) and the 
temperature structures (bottom panel) of the new (solid curves) 
and the old (dashed curves) models with the same g m dlg. The 
models are computed for pure hydrogen as well as solar abun- 
dance. The models have different effective temperatures, there- 
fore the spectra are normalized to the maximum flux and plotted 
against the scaled photon energy. For clarity, the spectra of so- 
lar abundance models are shifted by a factor of three, and the 
relative temperatures are shifted by adding one. 



ergies E/kT^, one also sees that the new spectra are harder (see 
Fig. |1. 

3.2. Limb-darkening 

Knowledge of the angular distribution of the emergent radiation 
is important when only part of the star is visible (for exam- 
ple, is partially blocked by the accretion disk), or when there 
are inhomogeneities at the NS surface related, for example, 
to the varying gravitational acceleration due to the rapid rota- 
tion. Computation of the amplitude of reflection from the ac- 
cretion disk also requires that information. The limb-darkening 
law (i.e. angular dependence of the intensity) for the radi- 
ation emerging from the optically thick electron-scattering- 
dominated atmosphere is described by the H^(p) function 
(IChan drasekh ar & Breed Il947t IChandrasekhar1ll960l; ISobolevi 
119491 119631) . A" simple approximation to this function is I(p) « 
1 +2.Q6p. Our simulations show that the intensity closely (within 
a couple of per cents) follows the H^Hp) function around the 
peak of the spectrum, in the observed energy range (see Fig . [5j ■ 
At photon energies below 1 keV, the computed angular distribu- 
tion becomes more isotropic, because there free-free absorption 
dominates over electron scattering and the upper atmosphere 
layers are almost isothermal. The low-inclination intensity is 
above and the high-inclination intensity is below the electron- 
scattering limb-darkening law at energies above the peak, be- 
cause the temperature of the layer where the photons originate 
drops with the inclination. The iron absorption edge at ~9 keV 
significantly affects the angular distribution for the solar abun- 
dance models (Fig. [5] bottom panel). Above the edge, radiation 
becomes more directed along the surface normal. 

The angular distribution of the intensity also depends on the 
specific form of the RF used in the calculations. The approxi- 
mate angle-dependent gives results very close to the reference 
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I = 0.8 models are multiplied by a factor of three for clarity. The 
lower sub-panels show the ratios of the corresponding models. 



intensity spectra, within 2% for the largest angle 6 (see Fig. [6] 
top panel). The exact, but angle-averaged, RF, does not give such 
accurate results (see Fig. O bottom panel). 

3.3. Color correction factors 

All computed emergent spectra of the new atmosphere models 
were fitted by a diluted blackbody spectrum 



F E * wB E (f c T eS ) 



(29) 



using five different fitting procedures described in SPW11. We 
calculated the color correction / c and dilution w factors in the 
energy band (3-20)x(l + z) keV corresponding to the observed 
range of the RXTEfPCA detector. We calculated redshifts from 
logg by adopting a NS mass equal to 1.4M (see Eqs. (Q]) and 
©): for logg = 14.0, 14.3, and 14.6, we getR = 14.80, 10.88, 
8.16 km and z = 0.18, 0.27, 0.42, respectively. Varying the mass 
in the interval 1-2M has a smaller than 0.1% effect on the color 
corrections (SPW11). The results of the fitting procedures are 
presented in Table[T](see also AppendixiDl). 

Deviations of the new model spectra from the best-fit diluted 
blackbodies are similar to those for the old spectra based on the 
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Fig. 6. Comparison of the emergent specific intensity at three 
angles (from the Gaussian quadrature p = 0.113,0.5,0.887, 
i.e. 6 = 83°5,60°,27!5) for the fiducial model computed us- 
ing the exact angle-dependent RF (solid curves in both panels) 
with the spectra based on the approximate angle-dependent RF 
(top panel, dashed curves) and exact angle-averaged RF (bottom 
panel, dashed curves). In the bottom sub-panels, the ratios of the 
accurate spectra to the approximate spectra for the three angles 
are shown. 



Kompaneets operator (Fig. |7). Comparison of the new and old 
color correction factors is shown in Figs. [8] and [9] Although the 
old f c — I dependences for different gravities were almost identical 
at high luminosities, we see that now these dependences deviate, 
because of the dependence of the limiting relative luminosity on 
\ogg. However, the old and new color-correction factors have 
very similar dependences on g ra d/g for all \ogg and chemical 
compositions (see Figs. [8]and|9]l, with the new / c being about 
1% larger. The d i fferen ce grows at g m dlg close to unity. 

iPavlov etaTI (1 19911) derived an approximate analytical for- 
mula for the ratio of the surface model atmosphere temperature 
to the effective one 



0.141n 3+5X +0.59 



-4/5 /3 + 5xX 2/15 



1 -I 



,3/20 



(30) 



which is correct for the highly luminous (/ > 0.9) model atmo- 
spheres. The numerical constants 0.14 and 0.59 were found from 
the fitting of their model atmospheres based on the Kompaneets 
operator description of Compton scattering. If the models were 
instead computed using an exact RF for Compton scattering, it 
is obvious that / should be substituted by the relative radiation 
acceleration g m dlg- We fit our color correction factors computed 
using the first fitting procedure (see SPW1 1) with a formula sim- 
ilar to Eq. ( 1301 and found numerical constants that also depend 
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Fig. 7. Relative deviations of the new, exact RF-based (solid 
curves) and old Kompaneets-based (dashed curves) spectra from 
the best-fit diluted blackbodies versus photon energy for hydro- 
gen (top panel) and solar H/He mixture with Z = 0.3Z o (bot- 
tom panel) low gravity (logg =14.0) models. Corresponding rel- 
ative luminosities and the effective temperatures are given at the 
curves. The vertical dotted line shows the lower boundary of the 
energy band, where the fitting procedures were performed. For 
clarity, models with I — 0.1 and 0.5 are shifted up by 0.2 and 0.1, 
respectively. 



on the chemical composition 



/ c * [0.102 + 0.008X] In 



3 + 5X 



1 -gmd/g 



+ 0.63 - 0.06X 



-4/5 



3 +5X 



2/15 



(grad/g) 



3/20 



(31) 



\ 1 - grad/g 

This approximation works well forg rac |/g > 0.8 (see Fig. ITOt. 
3.4. Radiative acceleration 

The radiative acceleration can formally be represented as a prod- 
uct of the flux and the temperature-dependent effective opacity 



°"s B r 4 f 



eff 



grad = K(T) 

c 

This expression can alternatively be written as 

grad _ l K(T) 
8 K e 



(32) 



(33) 
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Fig. 8. Color correction factors / c j (computed using method 1 
from SPW1 1) for model atmospheres of two chemical composi- 
tions (pure hydrogen and solar hydrogen/helium mix with 30% 
of solar heavy-element abundance) and different log g as func- 
tions of the relative luminosity I (top panel) or g ra d/g (bottom 
panel). The new models based on the exact RF are shown by 
the solid curves, and the old models, based on the Kompaneets 
operator, by dashed curves. 



In diffusion approximation, k(T) is given by the Rosseland mean 
opacity. When electron scattering dominates, it is often approxi - 
mated (neglecting the electron degeneracy) as (Paczvnski 1983) 



k(T) = k r (T) 



1 



kT 



38.8 keV 



0.86 



(34) 



This a pproximation is based on calculations bv lBuchler & Yuehl 
(119761) . who underestimated the opacity at low temperatures, 
where certain approximations were made because of severe nu- 
merical problems. A better approximation in the range 2-50 keV 
that is of interest here is (J. Poutanen et al., in preparation) 



k r (T) « K e 



1 + 



kT 



39.4 keV 



0.976 



(35) 



It is clear that the radiative acceleration decreases in the deep, 
hotter atmosphere layers and the ratio of radiation pressure force 
to the surface gravity decreases inwards (see Fig. ITTb . The ra- 
diative acceleration is smaller than that corresponding to the 
Thomson opacity, even in the upper atmosphere layers. The ac- 
tual radiative acceleration in the surface layers (see top panel 
in Fig. [TTI) computed from the models using Eq. d20l > is per- 
fectly described by Eqs. ( 1351 ) and ( 1331 throughout the whole at- 
mosphere (compare dotted and solid curves at the top panel in 
Fig. ITTl. while Paczynski's approximation underestimates it. 
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Table 1. Color-correction and dilution factors from the blackbody fits to the spectra of hydrogen atmosphere models at log g = 14.0. 



Set 



X=l Z = logg=14.0 T EM = 1.644 keV R = 14.80 km z = 0.18 



/ 


grad/g 


r cff (keV) 


/c,l 


/c,2 


/c,3 




fc,4 




fc,5 




W 2 / c 4 2 


"< 3 


0.001 


0.001 


0.292 


2.015 


2.009 


2.018 


1 


806 


1 


571 


0.698 


0.704 


0.695 


0.003 


0.003 


0.385 


1.830 


1.826 


1.833 


1 


662 


1 


473 


0.796 


0.800 


0.794 


0.010 


0.009 


0.520 


1.686 


1.682 


1.687 


1 


553 


1 


382 


0.877 


0.878 


0.876 


0.030 


0.028 


0.684 


1.589 


1.587 


1.589 


1 


525 


1 


375 


0.926 


0.926 


0.926 


0.050 


0.048 


0.777 


1.550 


1.550 


1.550 


1 


531 


1 


407 


0.942 


0.942 


0.942 


0.070 


0.067 


0.845 


1.531 


1.530 


1.530 


1 


533 


1 


431 


0.949 


0.949 


0.949 


0.100 


0.096 


0.924 


1.512 


1.511 


1.511 


1 


528 


1 


451 


0.954 


0.954 


0.954 


0.150 


0.143 


1.023 


1.493 


1.492 


1.492 


1 


512 


1 


462 


0.963 


0.962 


0.963 


0.200 


0.190 


1.099 


1.481 


1.479 


1.481 


1 


505 


1 


464 


0.960 


0.959 


0.960 


0.300 


0.285 


1.217 


1.481 


1.476 


1.480 


1 


500 


1 


478 


0.971 


0.969 


0.971 


0.400 


0.380 


1.307 


1.491 


1.486 


1.491 


1 


511 


1 


500 


0.972 


0.970 


0.972 


0.500 


0.474 


1.382 


1.503 


1.497 


1.503 


1 


521 


1 


515 


0.974 


0.972 


0.974 


0.550 


0.521 


1.416 


1.512 


1.506 


1.512 


1 


530 


1 


526 


0.975 


0.972 


0.975 


0.600 


0.567 


1.447 


1.522 


1.516 


1.522 


1 


539 


1 


538 


0.977 


0.973 


0.976 


0.650 


0.614 


1.476 


1.535 


1.530 


1.535 


1 


551 


1 


553 


0.979 


0.976 


0.979 


0.700 


0.661 


1.504 


1.551 


1.546 


1.551 


1 


566 


1 


569 


0.981 


0.977 


0.980 


0.730 


0.707 


1.530 


1.568 


1.563 


1.567 


1 


583 


1 


586 


0.982 


0.979 


0.982 


0.800 


0.733 


1.555 


1.589 


1.586 


1.588 


1 


600 


1 


608 


0.985 


0.983 


0.984 


0.850 


0.799 


1.578 


1.615 


1.613 


1.613 


1 


624 


1 


634 


0.987 


0.986 


0.987 


0.900 


0.844 


1.601 


1.646 


1.648 


1.644 


1 


653 


1 


665 


0.992 


0.993 


0.991 


0.950 


0.890 


1.623 


1.687 


1.694 


1.682 


1 


690 


1 


705 


0.997 


1.001 


0.996 


0.980 


0.917 


1.635 


1.718 


1.729 


1.712 


1 


718 


1 


735 


1.000 


1.007 


0.999 


1.000 


0.934 


1.644 


1.741 


1.756 


1.733 


1 


739 


1 


757 


1.002 


1.011 


1.000 


1.020 


0.952 


1.652 


1.772 


1.793 


1.762 


1 


767 


1 


787 


1.006 


1.018 


1.003 


1.040 


0.969 


1.660 


1.815 


1.844 


1.803 


1 


807 


1 


828 


1.010 


1.026 


1.007 


1.060 


0.981 


1.668 


1.880 


1.921 


1.864 


1 


869 


1 


888 


1.012 


1.036 


1.009 



Notes. Results for other chemical compositions and gravities are given in Table D. 1 available in electronic form at the Centre de Donnees as- 
tronomiques de Strasbourg (CDS). 




./) 



Fig. 9. Same as Fig. [8] but for different chemical compositions 
and logg = 14.0. 




Fig. 10. Color correction factor f c versus g m d/g for models of 
various chemical compositions. Solid curves are the results of 
calculations based on the exact Compton RF, the dotted curves 
give the approximation (fJTJ. For clarity, the curves for pure H 
models (shown for three logg) are shifted up by +0.2. 



bottom panel in Fig. ITTli. because the high-energy photons scat- 
ter on the relatively cold electrons in a predominantly forward 
direction, reducing the magnitude of the momentum transfer in 
comparison with the isotropic case. 

Radiative acceleration in the surface layers that we obtain 
from the models can be expressed through T e ff and / c as 



Radiative acceleration computed using the angle-averaged 
RFs is larger than that computed using angle-dependent RFs (see 



K (T = f c T eff ) = K e 



1 + 



kT 



38.8 keV 



(36) 
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Fig. 11. Top panel: Comparison of the relative radiative accel- 
eration in pure hydrogen atmospheres (solid curves) with our 
approximation given by Eq. ([35]) (dotted curves, which nearly 
coincide with solid curves) and with Paczynski's approximation 
(l34l > (dashed curves). Bottom panel: Comparison of the relative 
radiative acceleration for the fiducial model computed with dif- 
ferent RFs: exact and approximate, angle-dependent, and angle- 
averaged. The relative radiative acceleration for Thomson scat- 
tering is also shown. 



where a g = 1.01 +0.067(logg- 14.0) (see Fig.[T2>. The relation 
grad/g-l also slightly depends on the chemical composition, but 
the dependence on the surface gravity is stronger. 

The relation of the effective temperature to the effective 
opacity given by Eq. (f36l > allows us to estimate g m dlg for the 
given model parameters (I, T e ff) without actually computing the 
atmosphere models. We can make a first guess of / c , substitute 
T = f c T e s to Eq. ([36}, then use this k in Eq. d33l . find a new f c 
through Eq. (fJTJ), and iterate. 



4. Application to observations 

In our previous works (Sll and SPW11), we suggested a new 
method for the determination of NS radii and masses. It is based 
on the spectral (blackbody) normalization K at late phases of 
PRE X-ray bursts depending on the color correction factor only 



K s p bb (km) j 2 



D 



10 



1 

1 



fl(km)(l + z) 



(Af c T 



(37) 



where Djo = D/lOkpc is the distance. The observed relation 
K~ l ^—F\,\, between the blackbody flux Fbb and normalization K 
can be fitted by the theoretical dependence f c -l obtained from 
the atmosphere models. The fit gives two parameters: the value 
of A = [R (km) (1 + z)/D l0 ]~ 1/2 and the observed Eddington flux 
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Fig. 12. Dependence of the relative radiative acceleration g m &lg 
on I. Symbols represent the results of calculations based on the 
exact Compton RF, the solid curves give the approximation (f36t . 
The dot-dashed curve is the I - g m dlg relation. 



^Edd = ^Edd(l + zY 2 /(AnD 2 ). The distance-dependent quanti- 
ties A and F E dd can be combined to the distance-independent 
Eddington temperature which is the apparent effective tempera- 
ture corresponding to L E dd 



■ Edd.oo 



1/4 



1 



1 + Z 



6.4xl0 9 A^ d 4 d K. 



(38) 



If the distance to the source is known (for example, for sources 
in globular clusters), we can plot three curves on the M -R plane 
corresponding to the values of A, FEdd, and 7Edd,oo- Two crossing 
points give two pairs of NS mass and radius values that satisfy 
the observed data. 

It is now important to understand how the new models affect 
results based on formulae that use the Thomson cross-section 
for electron scattering. We have seen that owing to the Klein- 
Nishina reduction in the cross-section, the actual Eddington 
limit is reached at luminosities of 6-10% higher than LEdd 
and the shape of the f c —l relation differs somewhat from the 
Kompaneets-based results. As an illustration, we consider a long 
PRE X-ray burst of 4U 1724-307 in the globular cluster Terzan 2 
studied by Sll, who obtained a rather large NS radius for this 
source (R > 14 km). We now use new sets of theoretical rela- 
tions f a -l that we fit to the observed relation K ^—F^, taking 
pure hydrogen models as an example. We note that the new rela- 
tions significantly depend on the surface gravity (see Fig. [8] top 
panel), therefore, in principle, it would be possible also to find 
the surface gravity that provides the best fit to the observed rela- 
tion. This is, however, difficult in practise, because a change in 
log g can be compensated for by varying ^Edd- 

The theoretical f c -l curve for logg = 14.0 gives the best fit 
to the data (Fig. O for A = 0.168 and F E dd = 4.93 x 10~ 8 
erg s _1 cm 2 . The corresponding old, Kompaneets-based values 
are 0.170 and 5.25 x 10~ 8 , respectively. This small change in 
the best-fit parameters leads to a small decrease of TEdd.oo from 
1.64 x 10 7 K to 1.60 x 10 7 K and a corresponding increase in 
the NS radii (see Fig. [13] bottom panel). Interestingly, the new 
^Edd.oo curve crosses with the curve \ogg = 14.0 at the NS mass 
close to 1.5M , and the curve A = 0.168 also passes through the 
same crossing point when D\q = 0.53 (i.e. 5.3 kpc, the distance 
to Terzan 2. lOrtolani etd] [l997). The uncertainties in the ob- 
tained values arising owing to uncertainties in the data are very 
close to the differences between the old and the new values of 
the best-fit parameters (see Table 1 in Sll). 
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Fig. 13. Top panel: The fit of the X-ray burst data for 4U 1724- 
307 by the theoretical models for the NS atmosphere. The circles 
indicate the observed dependence of K~ 1 ' versus F for the long 
burst (S 1 1 ) and the solid curve corresponds to the best-fit theo- 
retical model for a pure hydrogen atmosphere and \ogg = 14.0. 
Vertical dotted and dashed lines show the best-fit value for F Edd 
for the old and new models. Bottom panel: Constraints on mass 
and radius of the NS in 4U 1724-307. The new (solid) and the 
old (dashed) curves corresponding to TEdd.oo are shown. The 
dash-dotted curve corresponds to the new best-fit parameter A 
for the distance to the source of 5.3 kpc. 



A few attempts to find NS mass and radius from the di- 
rect fitting of the observed X-ray burst s pectrum by NS at- 
mosphere model spectra were performed (|Maiczvn a & Madei 
120051: iMiller et al.ll201 ltlKusmierek et al.ll201 ll) . Thev have pro- 
posed finding the best fit for a fixed log g by varying the relative 
luminosity I (or effective temperature T s g ) and the gravitational 
redshift z. The best-fit between all trial log g values then gives 
the desired log g and z. Thereafter, the NS mass and radius can 
be found using Eqs. ([1) and (0). Unfortunately, the curves on 
the M-R plane corresponding to the fixed values of log g and 
z cross at very small angles (see Fig. 1 in SPW11). Therefore, 
the uncertainties are expected to be large. More importantly, the 
model spectra are very close to the blackbody in the observed 
RXTE/PCA energy band. This means that for an arbitrary prod- 
uct / c r e ff (i.e. color temperature of the theoretical spectrum) we 
can find the redshift z that makes the combination / c 7^/(1 + z) 
equal to the observed color temperature T^,, and in turn this 
method extremely unreliable. 



5. Summary 

We have considered hot NS model atmospheres taking into ac- 
count Compton scattering using the relativistic kinetic equa- 
tion with the exact angle-dependent RF and cross-section 



(Nagi rner & Poutanenl Il994t iPoutanen & Svenssonl Il996h . ac- 
counting also for the induced scattering. We have developed a 
method to solve the obtained radiation transfer equation using 
a short characteristic method with the accelerated A-iteration. 
We have implemented this solution in our computer code for the 
model atmosphere calculations (SPW1 1). 

We have examined the properties of the new model atmo- 
spheres. The main difference in comparison with the old models 
computed with the Kompaneets operator concerns radiative ac- 
celeration. In the new approach, the Klein-Nishina reduction in 
the electron scattering cross-section leads to a decrease in the 
radiative acceleration relative to that for the Thomson scatter- 
ing cross-section. The grid of model atmospheres was extended 
to higher effective temperatures up to luminosities formally ex- 
ceeding the Eddington limit Lndd computed for the Thomson 
cross-section. The role of the radiation pressure reduces in the 
deeper, hot, optically thick layers due to the same effect. 

We computed a new set of 484 hot NS model atmospheres. 
Following SPW11, we computed the models for six differ- 
ent chemical compositions (pure hydrogen, pure helium, solar 
hydrogen/helium mix with various abundances of heavy ele- 
ments Z = 1, 0.3, 0.1, and 0.01 Z Q ) and three surface gravities 
(logg = 14.0, 14.3, and 14.6). The relative luminosities range 
from /=0.001 to 1.06—1.1 (depending on \ogg). The new mod- 
els with I < 0.8 are almost identical to the old models based on 
the Kompaneets operator approach. At higher I, the deviations 
are more significant due to a different treatment of the radiative 
acceleration. The difference between the new and old models is 
small if instead of I one uses the relative radiative acceleration 
gmdlg as a parameter. 

The spectra of the new and old models deviate by less than 
5% in the region around the spectral peak, and the deviation 
grows at higher energies. We have also tested various approxi- 
mate and angle-averaged RFs for Compton scattering and found 
that they produce spectra that deviate from the exact solution 
typically by less than 2%. The compa rison of our spectra to the 
mode l s computed using Mad ei' s code ([Madei 1991: Madei et al.l 
l2004t iMajczvna etai1l2005l) revealed dramatic differences in 
the spectral shape. We attri bute these diffe rences to their us- 
age of the incorrect RF from Guilbert (119811) and/or to the non- 
convergence of their radiation transfer calculations. 

The spectra of all our models were fitted by the diluted black- 
body spectra in the RXTE/PCA energy band (3-20 keV) using 
the same five fitting procedures as in SPW1 1. The new f c -l rela- 
tions are similar in shape to the old relations, but depend signif- 
icantly on the surface gravity. However, the dependences of the 
color correction on g m dlg are almost indistinguishable, with the 
new / c being larger by approximately 1%. The color corrections 
with the corresponding dilution factors, the theoretical emergent 
spectral energy distributions, and specific intensities at various 
angles for all models are available at CDS. 

In the present paper, we have also tested how the new f c — 
I relations affect the determination of the NS masses and radii 
derived from the spectral evolution observed during the cool- 
ing stages of PRE bursts. We found that the best-fit NS radius 
increases by about 10% from the value obtained using the old 
model based on the Kompaneets equation (Sll). We note that 
bursting NSs are rapidly rotating and accounting for the latitu- 
dinal variation in gravity and the Doppler effect will affect the 
estimated NS radius. This will be a subject of a separate publi- 
cation. 
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Appendix A: Relativistic kinetic equation for 
Compton scattering and the RFs 

For completeness let us rederive here the exact relativis- 
tic expressions for the Compton scattering RFs. Fo r the 
detailed derivation, see iNagirner & Poutanenl d 19931) and 
iPoutanen & Vurml d2010l) . where more general problems have 
been solved. In the first paper, the redistribution matrices de- 
scribing Compton scattering of polarized radiation in terms of 
Stokes parameters were derived, while in the second the RFs 
for anisotropic electron distribution have been obtained. We start 
from the relativistic kinetic equation (RKE) for photons that de- 
scribes Compton scattering. 

A.1. Radiative transfer equation 

A description of interactions between photons and electrons via 
Compton scattering accounting for the induced scattering and 
electron degenerac y can be provided by the explicitly covariant 
RKE for photons (feuchler & Yueh| [l976t Ide Groot etaLlll980b 
INagirner & Poutanenll 19931 1 1994 



r 2 C dpdp x dx x . 
x ■ Vn(x) = — -r \ F Sr(p + x - p - x) 

24J r ri - 1 - 

x {n{xi)[l + n(X)]«e(/>i)[l - n e (p)} 
- n(x)[l + n(xi)]« e (»[l - n e (/>i)])} , 



(A.l) 



where V = {d/cdt, V} is the four-gradient, r e is the classical elec- 
tron radius, Ac = h/m e c is the Compton waveleng th, F is the 
Klein-Nishina reaction rate (Berestetskii et al. 1982) 



U ft/ U ft 

and 



ft ? 



^ = £1 ' -1 = E ' -' ft = Ei ' - = E ' -1 



(A.2) 



(A.3) 



are the four-products of the corresponding momenta (second 
equalities in Equations ( IA.3I ) arising from the four-momentum 
conservation law represented by the delta-function in Eq. ( IA.U ). 
Here we define the dimensionless photon four-momentum as 
x = {x,x} = x{l,<u}, where co is the unit vector in the photon 
propagation direction and x = hv/m e c 2 . The photon distribution 
is described by either the occupation number n or the specific 
intensity (per dimensionless energy interval) I{x) = x 3 n{x)/C, 
where the constant C is given by Eq. (1131 1. The dimensionless 
electron four- momentum is p = {y,p} - {y, pfl] = y{l,/3£l], 
where Q is the unit vector along the electron momentum, y and 
p = y/y 2 - 1 are the electron Lorentz factor and its momentum 
in units of m e c, and /3 is the velocity in units of c. The electron 
distribution is described by the occupation number n e . For the 
isotropic electron distribution, we use the electron distribution 
function f e (p) = 2n e (p)/ A c N e , normalized to unity 



■r 

Jo 



4n Up)p 2 dp=\. 



(A.4) 



In the following, we consider a steady state and ignore elec- 
tron degeneracy, because in the upper atmosphere layers, where 
the radiation spectrum is formed, electrons are non-degenerate. 
We define the RF as 



ion J y y\ - 1 



+ x l - p - x). (A.5) 
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For the rel afivistic Maxwellian dist ribution of temperature © 
kTJm e c 2 dJiittnerll 1 9 1 It ISvngdl 1 957h . 



1986; iPoutanen & Vurmll2010l) . choosing the polar axis along 
the direction of the transferred momentum 



Mp) = 



l 



4tt© K 2 (l/&) 



exp(-y/0) 



(A.6) 



(where K 2 is the modified Bessel function), the RF satisfies the 
symmetry property 



n = (x\(D\ - xa>) IQ, 
where 

Q 2 - Oi<Ui - xcj) 2 - (x - x\) 2 + 2q. 



(A. 14) 



(A. 15) 



R(x xi) e- x/@ =R(x 1 ^> x) e 



-Xi/S 



(A.7) 



Thus, the integration variables become cos a = n ■ £1 and the 
azimuth O. The RF ( I A. lit can then be written as 



which follows from its definition in Eq. iA.51 and the energy 
conservation y\ - y + x — x\, or from th e detailed balance 
condition (see eq. 8.2 in |Pomraninglll973t) . Using this result 
it is easy to show that the Bose-Einstein distribution n(x) = 
l/(exp{[x - yu]/©} - 1) with any chemical potential is a solution 

oftheRKE(EQ- 

In the absence of strong magnetic field, the medium is 
isotropic, therefore the RF depends only on the photon energies 
and the scattering angle (where rj is its cosine), i.e. we can write 
R(x\ — > x) — R(x, x\, if). The kinetic equation ( IA.U can then be 
recast in a standard form of the radiative transfer equation 

co ■ Vn(jr) 



2n 1 



R(x, x\,rj) = — — 

-f. Mp) Pdyj d<t> J 6(F) F dcosa, (A.16) 
1 o -1 

where now 

F = y(x\ - x) - q - pQcosa. (A. 17) 

Integrating over cos a using the delta-function, we get 



R(x, X\, rf) 



f 

Jy 



MP) R(x,x u t], y)dy, 



(A.18) 



+ 



where the integration over the electron distribution can been 
done numerically. Here we introduce the RF for monoenergetic 
electrons 



i r r 2 

= -n(x)— I x\dx\ I d a>\ R{x\,x,jf)[l + n(x\)} 
x Jo J 

i r r , 

[1 + n(x)]— J x\dx\jdco\R{x,x\,rj)n{x\). (A. 8) 

For the plane-parallel atmosphere, this reduces to 

P ' — =n(x,p)— I x\dx\ I d/j] R(x],/ii;x,fj)[l + n(xi,fi\)] 

dr T x Jo J_j 

i r°° r 1 

- [1 + n(x,fi)]— I xidxi I dfii R(x,/j.;x\,/j.\)n(xi,iJi), (A. 9) 
x Jo J-i 

where dTj = -o- T N e ds = K e dm, fj. and /Ji are the cosines of the 
angle relative to the normal, 77 = + yjl -p 2 -yjl - p\ cos <p y > y t {x,X\,rj) - (x - x\ + Q y/l + 2/qj /2. 



1 1 r n 

R( x ,x u i],y) = -— Fd<t>. 
Q2n Jo 

We need to substitute 
y(x\ -x)-q 



cos a — 



pQ 



(A. 19) 



(A.20) 



into the expression for F. The condition | cos a\ < 1 provides the 
constraint 



(A.21) 



and 

R(x 



R(x,xu 





7) dtp 



(A. 10) 



Integrating over azimuth <1> in Eq. ( IA. 191 > gives the exact analyti- 
cal expr ession for the RF that is valid for any photon and electron 
energy ( Aharonian & Atovanlll981l : lNa girner & Poutanenl fl994; 



is the azimuth-integrated RF. Rewriting Eq. JA.91 > in terms of the 
intensity I(x,fi), we get the radiative transfer equation ( TTUb that 
accounts for electron scattering with the scattering opacity and 
the source function given by Eqs. ( TTZb and (fT4b . respectively. 

A.2. Redistribution functions 



Pout anen & Vurmll2010h . which we use in our calculations 

- 2 -2q 



2 q" 

R(x,x U T], y) = — + — 



<t 



in 



1 + I 



d- d + ' 



(A.22) 



where 



(y-*) 2 + 1+?? 



The expression ( IA.5t for the RF can be simplified by tak- 
ing the integral over p with the help of the three-dimensional 
delta-function and using the identity 5{y\ + x\ — y — x) = 

y$\X] ■p 1 - x .-(p 1 



■a 2 ± 



Q 2 )/2. 



(y + Xl f + ^ 



(A.23) 



We n ote that the RF jA. 19b satisfies the detailed balance condi- 
ill 



R(x,x u t]) = ^- f ^-Mp)F8{T), 
167T J y 



(A.ll) 



tion (iNagirner & Po utane nlll994h 
R(x\, x, 77, y) — R(x, x\, 77, y + x — x\) 



(A.24) 



where we have dropped the subscript 1 from the electron quan- 
tities and 



r = y(x\ — x) — p(x\0)\ — xto) ■ Q - q, 
q = x ■ X\ — xx\{\ - rj). 



(A. 12) 
(A.13) 



To integrate over angles in Eq. JA.llll, we follo w the recipe 
proposed bv lAharonian & Atovanl (Il98ll) (see also lPrasad et all 



The RF (IA.1U is related to the scattering kernel (8.13) of 
iPomranina (1 1 973b as R(x, xi,rj) = <x s (xi — > x, ri)x\lx. The form 
given by Eq. ( IA. 161 > is equivalent to eq. (A4) of Buc hier & Yuehl 
(Il976h . The derived RF for monoenerge t ic elec trons ( IA.221 > is 
equivalent to eq. (A5) oflBuchler & Yuehl (fT976b and eq. (14) of 
I Aharonian & Atovanl (1 1 98 ll) . 

A very simple approximate expression for the RF can be 
obtained by assuming that the scattering in the electron rest 
frame proceeds in the Thomson regime (i.e. coherent) and is 
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isotropic. This is e quivalent to substituting F by 4/3 in Eq. 
( IA.19|. We then get (lArutvunvan & Nikogosvanll980tlPoutanenl 
1994: Poutanen & Svenssonl 19961) 



R{x,xi,tj,y) = — . 



(A.25) 



Integrating it over the Maxwellian distribution dA.61 > gives 

1 e-?'> & 



R(x, X[,T]) 



&rfiZ 2 (l/0) 



(A.26) 



This approximate RF is also used in the calculations. We note 
that this RF also satisfies the detailed balance condition dA.7l ). 



Appendix B: Method for solving the radiation 
transfer equation 

The formal solution of the radiative transfer equation gives a re- 
lation between the outward I + (x,/i) and the inward F(x,fi) = 
I(x, -fi) intensities at some depth point i (on the optical depth 
grid rf, i = 1, ...,N) with the adjacent intensities 



If(x,v) = /, + +1 (x, yU )exp[-(r+ 1 -Tt)/A + 

S + (t, x, fx) exp[-(f - rf)/fj] dt/fi, 
lj{x,n) = £ 1 (*,/0exp[-(Tr-T^ 1 )/|/u|] + 

S-(t,x,fx)exp[-(Tj -t)/\fx] dt/\n\. 



£ 



(B.l) 



(B.2) 



£ 



The integrals can be replaced by the sums using the parabolic 
approximation 



It(x,fi) = T^jC*, jU)exp(-AT+ +1 ) + 

afSUix,^ +/3lSf(x,n) + yfSf +l (x,ij), 
I~(x,fi) = /rj(^,ju)exp(-ATp + 

ajSj_ x (x,n) +/3JSJ(x,/j) + y~Sj +1 (x,fx), 

where 

Arf = (rf-rt_0/H 

_ _ _ g 2.,--( A7 7 + i +2AT D e u 
= < * + Ar r( A T;+I+ Ar r ) ' 
PI = ([Ar7 + Ar7 +I ] eu - e2 .)/(AT7Ar7 +1 ), 

r, r = ( e 2,;- AT7e u)/[ A W A ^ 7 + i +ATp], 
< = (^i-AT^^/CATttAT^+Arn), 
tf = ([Arf +A< 1 ]< /+1 -4 ;+I )/(ArfA< I ), 



(B.3) 
(B.4) 

(B.5) 



7i = 
and 



4; +1 -( A ^ + 2AT+ l)e + j+1 



1 - exp(-ATf ), 
e% = Arf -et ti , 

4,i = ( A ^) 2 "2<r 

At the first depth point, the coefficients are 
at = 0, # = < 2 /At+, y\ = el 2 -p{, 



(B.6) 



(B.7) 



£r=60keV, r=2 
AT = 1 keV 




x = hv/ m c 

Fig. B.l. Comparison of the emergent spectra from a hot elec- 
tron slab back-irradiated by a blackbody computed using the 
compps code and the code presented here. 



and at the last point they are 



K 0JV 



// At a 



7n 



0. 



(B.8) 



We note that the inward and the outward opacities along the 
same ray are different (see Eq. (fT2l ). 

The formal solution for a given source function starts for the 
inward intensities from the outer boundary condition (the lack of 
incoming radiation at the surface) 



h (-*>/") =yS I 5 1 (jc, / u) + y 1 5 2 (x, J u) 



(B.9) 



up to the last depth point N. The intensities at the innermost 
depth point are found using the inner boundary condition, which 
is taken from the diffusion approximation 



I^(x,fi) = I N (x,fi) + 2 



Ar« 



(B.10) 



The full solution is found iteratively using an accelerated A- 
iteration. At the first iteration, the thermal part of the source 
function is taken. For the subsequent iteration n, the intensities 
obtained from the previous iteration n - 1 are used to compute 
the current source functions Sf'". Iterations are continued un- 
til the relative change becomes smaller than the predetermined 
accuracy 



max 



J'Hx) 



J"~\x) 



1 



< 10" 



(B.ll) 



where 7,(x) are the mean intensities. This solution method of the 
radiation transfer equation was tested for a rather optically thin 
(Thomson optical depth tt = 2) and hot (kT e = 60 keV) electron 
slab back-illuminated by soft blackbody photons of kTsB — 1 
keV. The solution for the emergent intensities obtained at five an- 
gles using our method were compared wit h the solution obtained 
with the Comptonization code compps dPoutanen & Svenssonl 
19961 see Fig. ED. 

In the optically thick case (tj » 1), which is typical of NS 
atmospheres, the convergence of the solution can be accelerated 
using the following procedure. The difference between the for- 
mal solution obtained in the current iteration / ± n FS and the so- 
lution at iteration n - 1 is increased by some factor 



lf' n (x,fi)-lf' n l (x,p) = 



v few 



sf n {x,n) A* ; (x,ju) 



(B.12) 
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Fig. B.2. The maximum relative change in the solution of the ra- 
diation transfer equation at the final temperature correction com- 
puted by different versions of the accelerated A-iteration. 



where 

sf' n (x,n) 



<rf n (x,n) 
crf n (x,fi)+ki(xy 



(B.13) 



spectrum (solid curve) is much more peaked and softer than that 
computed by Madej's code (dashed curve). We also note that our 
spectrum is very close to the diluted blackbody spectra, which 
cannot be said about Madej's spectrum. 

We can suggest two hypotheses to explain the discrepancies. 
First, there is a difference in the RF used by the two codes. We 
employed the exact relativistic RF (see Appendix IA.2I ). which 
had been e xtensively studied and tested against M onte-Carlo 
simulations (Stern et al. 1995; Zdziarski et al. 2000). Madej and 
collab orators dMadeil 1 1 99 lb iMadei et al.ll2004t IMaiczvna et all 
2005) used the (angle -averaged) RF derived by Guilbert (1981) 
(see his Eqs. (8) and (10)), which differs from the correct expres- 
sion (|AT9]> by an additional factor [1 -Q3-n)(n-Wi)]/(l -fJ-a>i), 
which appeared from an error in the Jacobian. Guilbert's RF also 
does not satisfy the detailed balance condition ( IA.24l > and there- 
fore the RF integrated over the electron distribution does not sat- 
isfy the condition ( IA.7t either. Because the energy transfer by 
Compton scattering scales as ff, the error on the order of fj ob- 
viously invalidates all the results obtained with this RF. 

Second, some of the discrepancies could be connected to a 
difference in the computation of the radiation transfer and the 
atmosphere modeling. We illustrate this in Fig. IC.ll (top right 
panel). Using our radiation transfer code, we computed the spec- 
trum using the temperature structure from Madej's calculations. 



and A* ,(x, fx) is the diagonal term of the approximate A-operator when calculations proce eded until the accuracy of 10 



(B.14) 



x 2 r ^ r 1 

Jo x\ J_i 



dfiiR(x,/x;x u ^i). 



The acceleration is not high (about 30 - 40 %) and the number 
of necessary iterations is still large (see Fig. IB. 21 ). However, in 
the process of the model atmosphere computation it is possible 
to use the source function from the previous temperature itera- 
tion as the starting approximation for the current source function 
(see details in Section[2]). In this case, the acceleration depends 
on the value of the temperature corrections Ar, . At the first few 
temperature iterations, when AT 1 , are large, the acceleration is in- 
significant, but at later iterations, when Ar, are relatively small, 
the accelerated A-iterations converge very quickly (Fig. IB.2l ). 

Appendix C: Comparison with Madej's code 

The only other attempt to compute NS atmospheres using 
an integral approach to Compton scattering going beyond the 
Komp aneets approx i mation was that o f J. Madej and collabo - 
rators dMadeil 1 199 lb IMadei et al] 12004b IMaiczvna et all 120051) . 
They us ed an angle-ave raged RF for Compton scattering de- 
rived by Guilberi| ( 19811) . It is important to compare our results 
with those obtained by Madej's code for the same input param- 
eters. We selected our fiducial model as a testbed and computed 
two models. In the first, we approximated electron scattering 
by coherent Thomson scattering. In the second, the exact fully 
relativistic RF for Compton scattering was used. Calculations 
for identical parameters were performed with Madej's code (J. 
Madej and A. Rozanska, private communication). A comparison 
between the results is shown in Fig. lC. II 

We see that the Thomson scattering models are very close 
to each other, in terms of both the spectra and the tempera- 
ture structures (Fig. IC.ll left panels). However, the models with 
Compton scattering differ substantially (Fig. IC.ll right panels). 
Temperatures in the upper layers with column densities less than 
10 3 g cirT 2 are lower in our model by up to 5-15%, and our 



was 

reached (dotted curve), the spectrum was much above our bench- 
mark spectrum, i.e. it had a much higher effective (and color) 
temperature. If calculations were stopped at a lower accuracy of 
10~ 2 (dash-dotted curve), the spectrum was found to be simi- 
lar to Madej's spectrum. This suggests that the radiation transfer 
iterations with Madej's code did not sufficiently converge. 

Madej et al . use the partial linearization method described in 
Mihalasl d 1 9781 p. 179) and modified to include Compton scat- 
tering. Their radiation transfer equation is rewritten to include 
the linearized part of the radiative equilibrium equation. The so- 
lution of this generalized radiation transfer equation satisfies si- 
multaneously the radiative equilibrium to the first order. Using 
the computed radiation field, the temperature correction for the 
current temperature structure of the atmosphere is found. This 
procedure is iterated until some convergence criterion is satis- 
fied, for example, the emergent bolometric flux is accurate to 
within 0.1 %. In this method, the accuracy of the current so- 
lution of the radiation transfer equation is on the order AT/T, 
which is about 10 2 — 10 3 at the last iteration. We have shown 
above that this internal accuracy is insufficient to obtain the ex- 
act solution of the radiation transfer equation when Compton 
scattering is taken into account. We note that using this method 
is not possible to solve the radiation transfer equation for a given 
atmosphere model or, for example, for a homogeneous isother- 
mal slab. This means that it cannot be checked independently of 
the atmosphere modeling. 



Appendix D: Atmosphere model spectra and 
color-corrections 

Table D.l gives the color-correction and dilution factors from 
the blackbody fits to 484 atmosphere model spectra (fluxes), 
which in their turn are given in Table D.2. Table D.3 contains 
the emergent specific intensities at three angles. Tables D.l, 
D.2 and D3 are only available in electronic form at the CDS 
via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via 
http://cdsweb.u-strasbg.fr/cgi-bin/qcat7J/A-)- A/. 
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Fig. C.l. Left panels: Comparison of emergent spectra and temperature structures of the models computed by our code (solid 
curves) and Madej's code (dashed curves), when electron scattering is approximated by coherent Thomson scattering. Right panels: 
Same as left, but when the exact RF was used to compute Compton scattering. The spectrum computed by us using the Madej's 
temperature structure is shown by the dotted curve and the spectrum computed with a relative accuracy of 10~ 2 is presented by 
dash-dots. 
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